## ---- meetup-frequency ----

library(sp)
library(sf)
library(dplyr)

load('replication_file_02_meetup_event_venue.RData')
load('replication_file_03_spatial_analysis.RData')
load('replication_file_04_geocode_venue_italy.RData')

comuni_2018_sf <- 
  st_as_sf(comuni_sp)
comuni_2018_sf_buff_5k <- 
  st_buffer(st_transform(comuni_2018_sf, 32632), dist = 5000)
meetup_event$time <- 
  as.Date(meetup_event$time)
meetup_event_it_sf <- 
  meetup_event
meetup_event_it_sf$lon <- 
  meetup_venue_it@data$lon[
    match(meetup_event_it_sf$venue, meetup_venue_it$venue_id)] 
meetup_event_it_sf$lat <- 
  meetup_venue_it@data$lat[
    match(meetup_event_it_sf$venue, meetup_venue_it$venue_id)]
meetup_event_it_sf <- 
  meetup_event_it_sf[!is.na(meetup_event_it_sf$lon),]
meetup_event_it_sf <- 
  st_as_sf(meetup_event_it_sf, coords = c("lon", "lat"), 
           crs = 4326)
meetup_event_it_sf <- 
  st_transform(meetup_event_it_sf, 32632)

week_seq <- 
  as.character(seq(from = min(as.Date(meetup_event_it_sf$time)),
                   to = max(as.Date(meetup_event_it_sf$time)),
                   by = "1 weeks"))

meetup_7days_comune_8w_win <- data.frame()
for (this_date in week_seq[9:length(week_seq)]) {
  print(this_date)
  this_date <- as.Date(this_date)
  this_meetup_event_sf <- 
    meetup_event_it_sf %>% 
    dplyr::filter(yes_rsvp_count > 0 &
                    time > this_date - 8 * 7 &
                    time <= this_date)
  
  res <- st_join(this_meetup_event_sf, 
                 comuni_2018_sf_buff_5k, left = F) %>%
    as.data.frame() %>%
    select(event_id, PRO_COM_T, time)
  res$week <- cut(res$time, 
                  breaks = seq(from = this_date - 8 * 7, 
                               to = this_date + 7, "1 week"))
  res_sum <- 
    res %>% dplyr::group_by(PRO_COM_T) %>% 
    dplyr::summarize(meetup_7days = length(unique(week)))
  res_sum <- 
    merge(res_sum, 
          data.frame(PRO_COM_T = comuni_2018_sf_buff_5k$PRO_COM_T,
                     COD_PROV = comuni_2018_sf_buff_5k$COD_PROV,
                     COD_REG = comuni_2018_sf_buff_5k$COD_REG,
                     macro_region = comuni_2018_sf_buff_5k$macro_region,
                     stringsAsFactors = F), 
          all.y = T,
          by = 'PRO_COM_T')
  res_sum$meetup_7days[is.na(res_sum$meetup_7days)] <- 0
  res_sum$date <- this_date
  meetup_7days_comune_8w_win <- rbind(meetup_7days_comune_8w_win, res_sum)
}


meetup_7days_comune_8w_win$north_south <- 
  ifelse(meetup_7days_comune_8w_win$macro_region %in% 
           c('North-west', 'North-east',
             'Centre'),
         'North',
         'South')

meetup_7days_comune_8w_win_sum <- 
  meetup_7days_comune_8w_win %>%
  filter(PRO_COM_T %in% sample_comuni) %>%
  group_by(date, north_south) %>%
  summarize(mean = mean(meetup_7days),
            lower = t.test(meetup_7days)$conf.int[1],
            upper = t.test(meetup_7days)$conf.int[2])

meetup_7days_comune_8w_win$pop2011 <- 
  comuni_sp$pop2011[match(meetup_7days_comune_8w_win$PRO_COM_T, comuni_sp$PRO_COM_T)]

meetup_7days_treated_pop <-
  meetup_7days_comune_8w_win %>%
  group_by(north_south, date) %>%
  summarize(pop_treated = sum(pop2011[meetup_7days>1])/sum(pop2011))

meetup_7days_treated_pop_sample <-
  meetup_7days_comune_8w_win %>%
  group_by(north_south, date) %>%
  summarize(pop_treated = sum(pop2011[meetup_7days>1])/sum(pop2011))

meetup_7days_coeff <- data.frame()
for (this_date in as.character(unique(meetup_7days_comune_8w_win$date))) {
  print(this_date)
  this_date <- as.Date(this_date)
  this_df <- subset(meetup_7days_comune_8w_win, date == this_date)
  this_df <- merge(this_df, comuni_sp@data, by = 'PRO_COM_T')
  this_df$meetup_7days <- this_df$meetup_7days>1
  if (this_date <= as.Date('2013-02-24')) {
    lmmod_north <- 
      lm(M5S_2013~meetup_7days+log(pop_density)+log(dist_10k_ppl)+
           unemployment+housewive_perc+partwf_perc+over65_perc+
           degree_perc+foreignpop_perc+foreignpop_africa_perc+
           M5S_2013_neigh_avg+income2016_procapite+ANPIV_perc+ref16_turnout,
         data = this_df[this_df$north==T,])
    lmmod_south <- 
      lm(M5S_2013~meetup_7days+log(pop_density)+log(dist_10k_ppl)+
           unemployment+housewive_perc+partwf_perc+over65_perc+
           degree_perc+foreignpop_perc+foreignpop_africa_perc+
           M5S_2013_neigh_avg+income2016_procapite+ANPIV_perc+ref16_turnout,
         data = this_df[this_df$north==F,])
  } else {
    lmmod_north <- 
      lm(M5S_2018~meetup_7days+log(pop_density)+log(dist_10k_ppl)+
           unemployment+housewive_perc+partwf_perc+over65_perc+
           degree_perc+foreignpop_perc+foreignpop_africa_perc+
           M5S_2018_neigh_avg+income2016_procapite+M5S_2013+ANPIV_perc+ref16_turnout,
         data = this_df[this_df$north==T,])
    lmmod_south <- 
      lm(M5S_2018~meetup_7days+log(pop_density)+log(dist_10k_ppl)+
           unemployment+housewive_perc+partwf_perc+over65_perc+degree_perc+
           foreignpop_perc+foreignpop_africa_perc+M5S_2018_neigh_avg+
           income2016_procapite+M5S_2013+ANPIV_perc+ref16_turnout,
         data = this_df[this_df$north==F,])
  }
  coeff_mat_north <- summary(lmmod_north)$coefficients
  coeff_mat_south <- summary(lmmod_south)$coefficients
  
  meetup_7days_coeff <- 
    rbind(meetup_7days_coeff,
          data.frame(date = this_date,
                     north_south = c('North','South'),
                     meetup_7days_estimate = 
                       c(coeff_mat_north['meetup_7daysTRUE',1], 
                         coeff_mat_south['meetup_7daysTRUE',1]),
                     meetup_7days_stderr = 
                       c(coeff_mat_north['meetup_7daysTRUE',2], 
                         coeff_mat_south['meetup_7daysTRUE',2]),
                     meetup_7days_pvalue = 
                       c(coeff_mat_north['meetup_7daysTRUE',4], 
                         coeff_mat_south['meetup_7daysTRUE',4]),
                     stringsAsFactors = F))
}
meetup_7days_coeff$lower <- 
  meetup_7days_coeff$meetup_7days_estimate - 
  meetup_7days_coeff$meetup_7days_stderr
meetup_7days_coeff$upper <- 
  meetup_7days_coeff$meetup_7days_estimate + 
  meetup_7days_coeff$meetup_7days_stderr

save(meetup_7days_comune_8w_win, 
     meetup_7days_coeff, 
     file = "replication_file_11_meetup_7days_comune_8w_win.RData")


